home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #2 / Ham Radio 2000 - Volume 2.iso / HAMV2 / MISC / COIL200 / COIL.C next >
C/C++ Source or Header  |  1996-09-29  |  16KB  |  630 lines
  1. /*
  2.   Main program for electrical self inductance of a coil.
  3.  */
  4.  
  5. #include <stdlib.h>
  6. #include <stdio.h>
  7. #include <conio.h>
  8. #include <stdarg.h>
  9. #include <string.h>
  10. #include <ctype.h>
  11. #include <math.h>
  12.  
  13. #define SING_LAYER_CIRCULAR          '1'
  14. #define MULT_LAYER_CIRCULAR          '2'
  15. #define N_TURN_CIRCULAR_LOOP          '3'
  16. #define CIRCULAR_CURRENT_SHEET          '4'
  17. #define STRAIGHT_ROUND_WIRE          '5'
  18. #define TORUS_CIRCULAR_WINDING          '6'
  19. #define TORUS_RECTANGULAR_WINDING     '7'
  20. #define SING_LAYER_SQUARE          '8'
  21. #define MULT_LAYER_SQUARE          '9'
  22. #define WIRE_GAUGE              'G'
  23. #define INDUCTANCE              'I'
  24. #define CAPACITY              'C'
  25. #define FREQUENCY              'F'
  26. #define IMPEDANCE              'Z'
  27. #define UNITS                  'U'
  28.  
  29. #define PI 3.14159265358979323846
  30.  
  31. double floor(), sqrt(), exp(), log(), pow();
  32. double skin(), helsol(), nagaoka(), rsol(), squaresol();
  33.  
  34. double    diameter  = 0.0;
  35. double    length      = 0.0;
  36. double    dbundle      = 0.0;
  37. double    outer      = 0.0;
  38. double    frequency = 0.0;
  39. double    mu      = 1.0;
  40. double    L      = 0.0;
  41. double    cap      = 0.0;
  42. double    freq      = 0.0;
  43. double    w      = 0.0;
  44. double    N      = 0.0;
  45. double    gauge      = 0.0;
  46. double    units      = 0.001;
  47. char    compute     = SING_LAYER_CIRCULAR;
  48. int    pass_getnum;
  49.  
  50. extern double MACHEP; /* = 2.2e-16; */
  51. extern double MAXNUM; /* = 1.7e38; */
  52.  
  53. extern    int mono;
  54.  
  55. #ifdef __WATCOMC__
  56.   int  __cdecl textinit( void );
  57.   int  __cdecl textout( char *s);
  58. #else
  59.   int  textinit( void );
  60.   int  textout( char *s);
  61. #endif
  62.  
  63. char str[40];
  64.  
  65. int print( char *format, ... )
  66. {
  67.   char line[160];
  68.   va_list arglist;
  69.  
  70.   va_start( arglist, format );
  71.   vsprintf( line, format, arglist);
  72.   va_end( arglist );
  73.   return( textout( line ) );
  74. }
  75.  
  76. char *eng ( double val )
  77. {
  78.   static char *str = "-000.000 m";
  79.   int p, e, adp, ac;
  80.   char un;
  81.  
  82.   if ( fabs( val ) < 1e-13 )
  83.   {
  84.     val = 0;
  85.     p = 0;
  86.   }
  87.   else
  88.     p = (int) floor( log10( fabs( val ) ) );
  89.   e = p % 3;
  90.   if( e < 0 )
  91.     e += 3;
  92.   p -= e;
  93.   val *= pow( 10, -p );
  94.   p /= 3;
  95.   switch( p )
  96.   {
  97.     case -4:
  98.       un = 'p';
  99.       break;
  100.     case -3:
  101.       un = 'n';
  102.       break;
  103.     case -2:
  104.       un = 'u';
  105.       break;
  106.     case -1:
  107.       un = 'm';
  108.       break;
  109.     case  0:
  110.       un = 0;
  111.       break;
  112.     case  1:
  113.       un = 'K';
  114.       break;
  115.     case  2:
  116.       un = 'M';
  117.       break;
  118.     case  3:
  119.       un = 'G';
  120.       break;
  121.     case  4:
  122.       un = 'T';
  123.       break;
  124.     default:
  125.       un = '*';
  126.       break;
  127.   }
  128.   if( val >=0 )
  129.     val += 5e-4;    /* round after 3 didits */
  130.   ac = sprintf( str, "%d", (int) val );
  131.               /* 3 digits of precision */
  132.   if ( ( adp = (int) ( 1e3 * fabs( val - floor(val) ) ) ) != 0 )
  133.     ac += sprintf( str + ac, ".%03d", adp );
  134.   sprintf( str + ac, " %c", un );
  135.   return( str );
  136. }
  137.  
  138. /* Subroutine to get value entered on keyboard
  139.  * Displays previous value and updates only if a new value is entered.
  140.  */
  141. void getnum( char *s, char *su, double *pd, double units, int zflag )
  142. {
  143.   double t;
  144.   char *p;
  145.   char c;
  146.  
  147.   if( pass_getnum ) return;
  148.  
  149. gloop:
  150.   p = str + sprintf( str, "%s%s" , eng( *pd ), su );
  151.   while (*--p == ' ')
  152.     *p = 0;
  153.   print( "%%ff%s %%f2(%s)%%ff ? ", s, str );
  154.  
  155.   p = str;
  156.   while( (c = getch()) != '\r' && c != '\x1B' )
  157.   {
  158.     putch( c );
  159.     if( c == '\b' )
  160.     {
  161.       putch(' ');
  162.       putch('\b');
  163.       *--p = ' ';
  164.     }
  165.     else
  166.       *p++ = c;
  167.   }
  168.   *p = 0;
  169.   print( "\n" );
  170.  
  171.   if( c == '\x1B' )
  172.   {
  173.     pass_getnum = 1;
  174.     return;
  175.   }
  176.  
  177.   if( (str[0] == 0) && (*pd <= 0.0) )
  178.   {
  179.     t = *pd;
  180.     goto tests;
  181.   }
  182.  
  183.   if( str[0] != 0 )
  184.   {
  185.     sscanf( str, "%lf", &t );
  186.   tests:
  187.     if( zflag >= 0 )
  188.     {
  189.       if( t < 0.0 )
  190.       {
  191.     print( "%%fc ? Negative value not allowed; try again.\n" );
  192.     goto gloop;
  193.       }
  194.       if( zflag )
  195.       {
  196.     if( t == 0.0 )
  197.     {
  198.       print( "%%fc ? Zero value not allowed; try again.\n" );
  199.       goto gloop;
  200.     }
  201.       }
  202.     }
  203.     *pd = t * units;
  204.   }
  205.   sprintf( str, "%s%s" , eng( *pd ), su );
  206.   p = str + strlen(str);
  207.   while (*--p == ' ')
  208.     *p = 0;
  209.   print( "%%fb%s\n", str);
  210. }
  211.  
  212. /* Skin effect correction.
  213.  * This is a close empirical approximation to a table
  214.  * found in the CRC handbook.
  215.  */
  216. double skin( double d )      /* wire diameter */
  217. {
  218.   double s;
  219.  
  220.   getnum( "Frequency for skin effect [MHz]", "Hz", &frequency, 1e6, 0 );
  221.   if( frequency <= 0.0 )
  222.     return( 0.25 );
  223.  
  224.   s = 0.1071 * d * sqrt(frequency);
  225.   if( s > 100.0 )
  226.     s = 7.088/s;
  227.   else
  228.   {
  229.     s = 0.00546 * pow(s,5.0);
  230.     s = 0.25 * pow( 1.0+s, -0.2 );
  231.   }
  232.   return(s);
  233. }
  234.  
  235. main(int argc,char *argv[])
  236. {
  237.   double b, p, d, l, a, s;
  238.   char c;
  239.   char *Z;
  240.  
  241.   if( !textinit( ) )
  242.   {
  243.     printf( "\nCouldn't initialise text mode" );
  244.     printf( "\nExiting ..." );
  245.     exit( 0 );
  246.   }
  247.   if( !mono )
  248.   {
  249.     if( argc > 1 )
  250.       if( *argv[1] == '/' || *argv[1] == '-' )
  251.       {
  252.     c = tolower( *(argv[1] + 1) );
  253.     if( c == 'b' )
  254.       mono = 1;
  255.       }
  256.   }
  257.   setbuf( stdout, NULL );
  258.  
  259.   print( "%%ae0\n INDUCTANCE CALCULATOR %%b0\n" );
  260.   print( "%%f4by Steve MOSHIER and Emil LAURENTIU (YO3GGH)\n" );
  261.   if( !mono )
  262.     print( "%%f7Use /b on command line for black and white\n" );
  263.  
  264.   print( "%%fb\nDimension unit set to millimeter\n" );
  265.   units = 0.001;
  266.  
  267. loop:
  268.   pass_getnum = 0;
  269.   print( "%%ff\nChoose operation:\n" );
  270.   print( "\n%%fe%c%%f2 %-45s %%fe%c%%f2   %s", SING_LAYER_CIRCULAR,
  271.     "Single-layer circular solenoid of round wire",
  272.     WIRE_GAUGE, "Wire gauge calculation" );
  273.   print( "\n%%fe%c%%f2 %-45s", MULT_LAYER_CIRCULAR,
  274.     "Multilayer circular solenoid" );
  275.   print( "\n%%fe%c%%f2 %-45s %%fe%c%%f2   %s", N_TURN_CIRCULAR_LOOP,
  276.     "N-turn circular loop",
  277.     INDUCTANCE, "Inductance at resonance" );
  278.   print( "\n%%fe%c%%f2 %-45s %%fe%c%%f2   %s", CIRCULAR_CURRENT_SHEET,
  279.     "Circular solenoidal current sheet",
  280.     CAPACITY, "Capacity at resonance" );
  281.   print( "\n%%fe%c%%f2 %-45s %%fe%c%%f2   %s", STRAIGHT_ROUND_WIRE,
  282.     "Straight round wire",
  283.     FREQUENCY, "Frequency at resonance" );
  284.   print( "\n%%fe%c%%f2 %-45s %%fe%c%%f2   %s", TORUS_CIRCULAR_WINDING,
  285.     "Circular toroid, circular winding",
  286.     IMPEDANCE, "Impedance of LC" );
  287.   print( "\n%%fe%c%%f2 %-45s", TORUS_RECTANGULAR_WINDING,
  288.     "Circular torus ring, rectangular winding" );
  289.   print( "\n%%fe%c%%f2 %-45s %%fe%c%%f2   %s", SING_LAYER_SQUARE,
  290.     "Single-layer square solenoid",
  291.     UNITS, "Change unit for dimensions" );
  292.   print( "\n%%fe%c%%f2 %-45s", MULT_LAYER_SQUARE,
  293.     "Multilayer square solenoid (low precision)" );
  294.   print( " %%feEsc%%f2 Exit" );
  295.   print( "\n? (%%fe%c%%f2) ", compute );
  296.  
  297.   c = toupper( getche( ) );
  298.   if( c != '\r' )
  299.     compute = c;
  300.   print( "\n\n" );
  301.   switch ( compute )
  302.   {
  303.     case '\x1B':
  304.       print( "%%f7Exiting. Have a nice day!\n" );
  305.       exit(0);     /* graceful exit to monitor */
  306.  
  307.     case UNITS:
  308.     getunits:
  309.       print( "%%ff\nChoose the unit for dimensions:\n" );
  310.       print( "%%fe\nM%%f2 Millimeter" );
  311.       print( "%%fe\nC%%f2 Centimeter" );
  312.       print( "%%fe\nI%%f2 Inch\n" );
  313.       c = toupper( getche( ) );
  314.       print( "\n\n" );
  315.       switch ( c )
  316.       {
  317.     case 'M':
  318.       print( "%%fbDimension unit set to millimeter\n" );
  319.       units = 0.001;
  320.       break;
  321.     case 'C':
  322.       print( "%%fbDimension unit set to centimeter\n" );
  323.       units = 0.01;
  324.       break;
  325.     case 'I':
  326.       print( "%%fbDimension unit set to inch\n" );
  327.       units = 0.0254;
  328.       break;
  329.     default:
  330.       print( "%%fcInvalid answer\n" );
  331.       goto getunits;
  332.       }
  333.       goto loop;
  334.  
  335.     case WIRE_GAUGE:
  336.       getnum( "American Wire Gauge number", "B&S", &gauge, 1, -1 );
  337.  
  338.       d = 0.32485 * exp( -0.11594 * gauge );
  339.  
  340.       print( "%%fcWire diameter = %%fb%.3e inches (%.3f mm)\n", d, 25.4 * d );
  341.       d *= 2.54; /* convert to centimeters */
  342.       getnum( "Coil length", "m", &length, units, 0 );
  343.       getnum( "Winding thickness", "m", &dbundle, units, 0 );
  344.       if( pass_getnum ) goto loop;
  345.       b = 100 * length;
  346.       c = 100 * dbundle;
  347.       N = 1.0;
  348.       b = floor( b/d );
  349.       if( b > 0.0 )
  350.     N = b;
  351.       c = floor( c/d );
  352.       if( c > 0.0 )
  353.     N *= c;
  354.       print( "%%fc%.0lf turns will fit in this space.\n", N );
  355.       goto loop;
  356.  
  357.     case CAPACITY:
  358.       getnum( "Inductance value [nH]", "H", &L, 1e-9, 1 );
  359.       getnum( "Resonance frequency [MHz]", "Hz", &freq, 1e6, 1 );
  360.       if( pass_getnum ) goto loop;
  361.       w = 2 * PI * freq;
  362.       cap = 1 / ( w * w * L );
  363.       print( "%%fc\nThe capacity = %%fb%sF\n", eng( cap ) );
  364.       goto loop;
  365.  
  366.     case INDUCTANCE:
  367.       getnum( "Capacitor value [pF]", "F", &cap, 1e-12, 1 );
  368.       getnum( "Resonance frequency [MHz]", "Hz", &freq, 1e6, 1 );
  369.       if( pass_getnum ) goto loop;
  370.       w = 2 * PI * freq;
  371.       L = 1 / ( w * w * cap );
  372.       print( "%%fc\nThe inductance = %%fb%sH\n", eng( L ) );
  373.       goto loop;
  374.  
  375.     case FREQUENCY:
  376.       getnum( "Inductance value [nH]", "H", &L, 1e-9, 1 );
  377.       getnum( "Capacitor value [pF]", "F", &cap, 1e-12, 1 );
  378.       if( pass_getnum ) goto loop;
  379.       freq = 1 / ( 2 * PI * sqrt ( L * cap ) );
  380.       print( "%%fc\nThe resonance frequency = %%fb%sHz\n", eng( freq ) );
  381.       goto loop;
  382.  
  383.     case IMPEDANCE:
  384.       getnum( "Inductance value [nH]", "H", &L, 1e-9, 1 );
  385.       getnum( "Capacitor value [pF]", "F", &cap, 1e-12, 1 );
  386.       getnum( "Resonance frequency [MHz]", "Hz", &freq, 1e6, 1 );
  387.       if( pass_getnum ) goto loop;
  388.       w = 2 * PI * freq;
  389.       Z = strdup( eng( w * L - 1 / (w * cap) ) );
  390.       print( "%%fc\nThe LC impedance (at %sHz) = %%fb%sohm\n", eng( freq ), Z );
  391.       free( Z );
  392.       goto loop;
  393.  
  394.     case MULT_LAYER_CIRCULAR:    /* Thick circular solenoid */
  395.       getnum( "Average diameter", "m", &diameter, units, 0 );
  396.     getrad:
  397.       getnum( "Radial depth of winding", "m", &dbundle, units, 0 );
  398.       if( dbundle > diameter )
  399.       {
  400.     print( "%%fcDepth cannot exceed diameter. Try again.\n" );
  401.     goto getrad;
  402.       }
  403.       getnum( "Length of winding", "m", &length, units, 0 );
  404.       getnum( "Total number of turns", "", &N, 1, 1 );
  405.       d = 100 * diameter;
  406.       l = 100 * length;
  407.       a = 100 * dbundle;
  408.       L = rsol( 0.5*d, l, a, N );
  409.       break;
  410.  
  411.     case TORUS_RECTANGULAR_WINDING:
  412.     /* Circular torus ring of rectangular cross section
  413.      * Ref: Skilling 1948
  414.      */
  415.       getnum( "Inner diameter of circular torus ring", "m",
  416.           &diameter, units, 1 );
  417.       getnum( "Radial width of winding cross section", "m",
  418.           &dbundle, units, 0 );
  419.       getnum( "Height of winding cross section", "m", &length, units, 0 );
  420.       getnum( "Total number of turns", "", &N, 1, 0 );
  421.       d = 100 * diameter;
  422.       l = 100 * length;
  423.       a = 100 * dbundle;
  424.       d *= 0.5;
  425.       L = 2.0e-9 * N * N * l * log( 1.0 + a/d );
  426.       break;
  427.  
  428.     case MULT_LAYER_SQUARE:
  429.       getnum( "Side of square, to center of winding thickness", "m",
  430.            &diameter, units, 1 );
  431.       getnum( "Winding length", "m", &length, units, 1 );
  432.       getnum( "Winding thickness", "m", &dbundle, units, 0 );
  433.       getnum( "Total number of turns", "", &N, 1, 0 );
  434.  
  435.       d = 100 * diameter;
  436.       l = 100 * length;
  437.       a = 100 * dbundle;
  438.  
  439.       if( (a+l) == 0.0 )
  440.     goto operror;
  441.       s = d/(a+l);
  442.       L = 8.e-9 * d * N * N * (log(s) + 0.2235/s + 0.726 );
  443.       break;
  444.  
  445.     case SING_LAYER_SQUARE:
  446.     /* Square coil, rectangular winding cross section (thickness, length).
  447.      * Ref: Grover 1922a, CRC handbook
  448.      */
  449.       getnum( "Side of square, to center of winding thickness", "m",
  450.            &diameter, units, 1 );
  451.       getnum( "Winding length", "m", &length, units, 1 );
  452.       getnum( "Total number of turns", "", &N, 1, 0 );
  453.  
  454.       d = 100 * diameter;
  455.       l = 100 * length;
  456.       a = 100 * dbundle;
  457.  
  458.       L = squaresol( d, l, N );
  459.       break;
  460.  
  461.     case TORUS_CIRCULAR_WINDING:
  462.     /* Circular toroid
  463.      * Ref: CRC handbook
  464.      */
  465.     torrad:
  466.       getnum( "Radius distance from center of torus to center of winding", "m",
  467.           &diameter, units, 0 );
  468.       getnum( "Diameter of circular winding", "m", &length, units, 0 );
  469.       if( length > 2.0*diameter )
  470.       {
  471.     print( "%%fcDiameter cannot exceed toroid diameter. Try again.\n" );
  472.     goto torrad;
  473.       }
  474.       getnum( "Total number of turns", "", &N, 1, 0 );
  475.  
  476.       d = 100 * diameter;
  477.       l = 100 * length;
  478.  
  479.       L = 4.0e-9 * PI * N * N * (d - sqrt(d*d - 0.25*l*l));
  480.       break;
  481.  
  482.     case CIRCULAR_CURRENT_SHEET:
  483.     case SING_LAYER_CIRCULAR:
  484.     case STRAIGHT_ROUND_WIRE:
  485.     case N_TURN_CIRCULAR_LOOP:
  486.     /* Solenoid, loop, or straight wire
  487.        */
  488.       if( compute == STRAIGHT_ROUND_WIRE )     /* straight wire */
  489.     diameter = 0.0;
  490.       else
  491.     getnum( "Diameter of winding", "m", &diameter, units, 0 );
  492.  
  493.       if( compute == N_TURN_CIRCULAR_LOOP )
  494.       { /* circular loop */
  495.     length = 0.0;
  496.       }
  497.       else
  498.       {
  499.     if( compute == STRAIGHT_ROUND_WIRE )
  500.       getnum( "Length of wire", "m", &length, units, 0 );
  501.     else
  502.       getnum( "Length of winding", "m", &length, units, 0 );
  503.       }
  504.  
  505.       if( (length == 0.0) && (diameter == 0.0) )
  506.       {
  507.     L = 0.0;
  508.     break;
  509.       }
  510.  
  511.       if( compute == SING_LAYER_CIRCULAR )
  512.     getnum( "Wire diameter", "m", &dbundle, units, 1 );
  513.  
  514.       d = 100 * diameter;
  515.       l = 100 * length;
  516.       a = 100 * dbundle;
  517.  
  518.       if( d == 0.0 )
  519.       {
  520.       /* straight */
  521.       /* Ref: CRC handbook */
  522.  
  523.     print( "%%faStraight piece of wire.\n" );
  524.     getnum( "Please enter diameter of wire", "m", &dbundle, units, 1 );
  525.  
  526.     a = 100 * dbundle;
  527.     print( "%%faStraight piece of wire.\n" );
  528.  
  529.     L = 2.0e-9 * l * ( log(4.0*l/a) - 1.0 + skin(a) + 0.5*a/l );
  530.     break;
  531.       }
  532.  
  533.       getnum( "Total number of turns", "", &N, 1, 0 );
  534.  
  535.       if( l == 0.0 )
  536.       {
  537.       /* circloop */
  538.       /* Ref: CRC handbook; Snow 1952 */
  539.  
  540.     print( "%%faCircular loop.\n" );
  541.     getnum( "Please enter diameter of wire or bundle of wires", "m",
  542.          &dbundle, units, 1 );
  543.  
  544.     a = 100 * dbundle;
  545.  
  546.     p = d/a;
  547.       /* For "natural" current distribution b = -1
  548.        * For uniform current distribution b = 0 (Litz wire)
  549.        */
  550.     b = -1.0;
  551.  
  552.     if( N <= 1.0 )
  553.       s = skin(a);
  554.     else
  555.       s = 0.25;
  556.  
  557.     if( a <= 0.0 )
  558.       L = MAXNUM;
  559.     else
  560.     {
  561.       L = 2.e-9 * PI * N * N * d
  562.         * ( (1.0 + (2.0*b+1.0)/(8.0*p)) * log(8.0*p) - 2.0 + s
  563.         + (b-1.0)*(b-2.0/3.0)/(16.0*p*p));
  564.     }
  565.     break;
  566.       }
  567.  
  568.  
  569.       if( compute == CIRCULAR_CURRENT_SHEET )
  570.     L = nagaoka( 0.5*d, l, N );
  571.       else
  572.       {
  573.     if( a*N > l + 1e-3 ) /* allow a 10 um tolerance */
  574.     {
  575.       print( "%%fcWire diameter too large.\n" );
  576.       goto operror;
  577.     }
  578.     L = helsol( 0.5*d, l, a, N );
  579.       }
  580.       break;
  581.  
  582.     default:
  583.       goto operror;
  584.   }
  585.   if( pass_getnum ) goto loop;
  586.  
  587.   print( "%%fc\nThe inductance = %%fb%sH\n", eng( L ) );
  588.  
  589. fcloop:
  590.  
  591.   print( "%%ff\nCompute resonance:\n" );
  592.   print( "%%fe\nEsc%%f2 None" );
  593.   print( "%%fe\nM%%f2   Set mu for core (default 1)" );
  594.   print( "%%fe\nF%%f2   Capacity  -> Frequency" );
  595.   print( "%%fe\nC%%f2   Frequency -> Capacity" );
  596.   print( "\n?   (%%feEsc%%f2) " );
  597.   c = toupper( getche( ) );
  598.   print( "\n\n" );
  599.   switch ( c )
  600.   {
  601.     case '\x1B':
  602.     case '\r':
  603.       goto loop;
  604.     case 'M':
  605.       getnum( "Set the equivalent mu for inductance", "", &mu, 1, 1 );
  606.       print( "%%fc\nThe equivalent inductance = %%fb%sH\n", eng( mu * L ) );
  607.       break;
  608.     case 'F':
  609.       getnum( "Capacitor value [pF]", "F", &cap, 1e-12, 1 );
  610.       freq = 1 / ( 2 * PI * sqrt( mu * L * cap ) );
  611.       print( "%%fc\nThe resonance frequency = %%fb%sHz\n", eng( freq ) );
  612.       break;
  613.     case 'C':
  614.       getnum( "Resonance frequency [MHz]", "Hz", &freq, 1e6, 1 );
  615.       w = 2 * PI * freq;
  616.       cap = 1 / ( w * w * mu * L );
  617.       print( "%%fc\nThe capacitor value = %%fb%sF\n", eng( cap ) );
  618.       break;
  619.     default:
  620.       print( "%%fcInvalid answer\n" );
  621.       goto fcloop;
  622.     }
  623.   goto fcloop;
  624.  
  625. operror:
  626.   print( "%%fc? operator error\n" );
  627.   goto loop;
  628.  
  629. }
  630.